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We characterize the structure of maximally random jammed (MRJ) sphere packings by computing 
the Minkowski functionals (volume, surface area, and integrated mean curvature) of their associ¬ 
ated Voronoi cells. The probability distribution functions of these functionals of Voronoi cells in 
MRJ sphere packings are qualitatively similar to those of an equilibrium hard-sphere liquid and 
partly even to the uncorrelated Poisson point process, implying that such local statistics are rela¬ 
tively structurally insensitive. This is not surprising because the Minkowski functionals of a single 
Voronoi cell incorporate only local information and are insensitive to global structural information. 
To improve upon this, we introduce descriptors that incorporate nonlocal information via the corre¬ 
lation functions of the Minkowski functionals of two cells at a given distance as well as certain cell-cell 
probability density functions. We evaluate these higher-order functions for our MRJ packings as 
well as equilibrium hard spheres and the Poisson point process. It is shown that these Minkowski 
correlation and density functions contain visibly more information than the corresponding standard 
pair-correlation functions. We find strong anticorrelations in the Voronoi volumes for the hyperuni¬ 
form MRJ packings, consistent with previous findings for other pair correlations [A. Donev et al, 
Phys. Rev. Lett. 95, 090604 (2005)], indicating that large-scale volume fluctuations are suppressed 
by accompanying large Voronoi cells with small cells, and vice versa. In contrast to the aforemen¬ 
tioned local Voronoi statistics, the correlation functions of the Voronoi cells qualitatively distinguish 
the structure of MRJ sphere packings (prototypical glasses) from that of not only the Poisson point 
process but also the correlated equilibrium hard-sphere liquids. Moreover, while we did not find 
any perfect icosahedra (the locally densest possible structure in which a central sphere contacts 12 
neighbors) in the MRJ packings, a preliminary Voronoi topology analysis indicates the presence of 
strongly distorted icosahedra. 


I. INTRODUCTION 

Packings of frictionless monodisperse hard spheres in 
three dimensions serve as a simple, yet effective tool for 
modeling the complex behavior of such diverse many- 
particle systems as crystals, colloids, liquids, glasses, het¬ 
erogeneous materials, foams, and biological systems [1- 
10]. Among the rich multitude of states they are known 
to exhibit, considerable interest has been given towards 
sphere packings that are jammed (i.e., mechanically sta¬ 
ble), including maximally dense packings, low-density 
crystals, and amorphous packings [11-20]. 

In order to characterize the properties of sphere pack¬ 
ings, one may employ a geometric-structure approach 
in which configurations are considered independently of 
both their frequency of occurrence and the algorithm by 
which they are created [19]. For example, the simplest 
characteristic of a sphere packing is its packing fraction 
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(p, i.e., the fraction of space occupied by the spheres. 
Other useful characteristics of the structure include its 
pair-correlation function [14, 21-28], the pore-size distri¬ 
bution [4, 29], and structure factor [3, 24, 30-32]. 

It is also valuable to quantify the degree of ordering 
in a packing, especially those that are jammed (defined 
more precisely below). To this end, a variety of scalar 
order metrics ip have been suggested [19, 33] in which 
Ip = 0 is defined as the most disordered state (i.e. the 
Poisson point process) and "(/i = 1 is the most ordered 
state. Using the geometric-structure approach, one may 
therefore construct an “order map” in the (p-ip plane [19], 
where the jammed packings form a subset in this map. 
The boundaries of the jammed region are optimal in some 
sense, including, for example, the densest packings (the 
face-centered-cubic crystal and its stacking variants with 
'/'max = tt/vTS « 0.74 [13]) and the least dense jammed 
packings (conjectured to be the “tunneled crystals” with 
(pmin = 2(/)max/3 [16]). 

Among the set of all isotropic and statistically homo¬ 
geneous jammed sphere packings, the maximally random 
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jammed (MRJ) state is that which minimizes some order 
metric This definition makes mathematically precise 
the familiar notion of random closed packing (RCP) in 
that it can be unambiguously identified for a particular 
choice of order metric. A variety of sensible, positively 
correlated order metrics produce an MRJ state with the 
same packing fraction 0.64 [34, 35], which agrees roughly 
with the commonly suggested packing fraction of RCP 
in three dimensions [19]. However, we stress that density 
alone is not sufficient to characterize a packing; in fact, 
packings with a large range of "0 may be observed at this 
packing fraction [27, 34]. 

In order to study the MRJ state, a precise definition 
of jamming is needed. Therefore, rigorous hierarchical 
jamming categories have been defined [36, 37]: A pack¬ 
ing is locally jammed if no particle can move while the 
positions of the other particles are fixed; it is collectively 
jammed if no subset of particles can move without de¬ 
forming the system boundary; and if also a deformation 
of the system boundary is not possible without increas¬ 
ing its volume, the packing is strictly jammed^ i.e., it is 
stable against both uniform compression and shear defor¬ 
mations [38]. Strictly jammed MRJ sphere packings often 
contain a small fraction of rattlers (unjammed particles), 
but the remainder of the packing, i.e., the mechanically 
rigid backbone, is strictly jammed [39]. 

Determining the contact network of a packing is a 
subtle problem requiring high numerical fidelity. The 
Torquato-Jiao (TJ) sphere packing algorithm [25] meets 
this challenge by efficiently producing highly disordered, 
strictly jammed packings with unsurpassed numerical fi¬ 
delity as well as ordered packings [25, 40]. The algo¬ 
rithm achieves this by solving a sequence of linear pro¬ 
grams which iteratively densify a collection of spheres 
in a deformable periodic cell subject to locally lin¬ 
earized nonoverlap constraints. The resulting packings 
are, by definition, strictly jammed and they are with 
high probability exactly isostatic (meaning that they pos¬ 
sess the minimum number of contacts required for jam¬ 
ming) [39, 41]. The TJ packing protocol is intrinsi¬ 
cally capable of producing MRJ states with very high 
fidelity [25, 39]. The MRJ state can be regarded as a 
prototypical glass because it is maximally disordered (ac¬ 
cording to a variety of order metrics) while having infi¬ 
nite elastic moduli [19]. Atkinson et al. [39] recently car¬ 
ried out a detailed characterization of the rattler popula¬ 
tion in these MRJ sphere packings. They found a rattler 
fraction of 1.5% (substantially lower than other packing 
protocols, such as the well-known Lubachevsky-Stillinger 
packing algorithm [42]) and that the rattlers were highly 
spatially correlated, implying that they have a significant 
influence on the structure of the packing [29]. Moreover, 
as in previous studies [41], it was shown [39] that the 
backbone of the MRJ state is isostatic. We include rat¬ 
tlers in our analysis unless stated otherwise. 

MRJ packings have been characterized using a vari¬ 
ety of statistical descriptors, including the radial pair- 
correlation function 52 (^) (p ^52 is the probability den- 



FIG. 1. (Color online) Maximally random jammed (MRJ) 
sphere packing and its Voronoi diagram. Among all jammed 
sphere packings (roughly speaking, the mechanically stable 
packings), the MRJ state is the most disordered one. 


sity for finding two sphere centers separated by a ra¬ 
dial distance r, where p is the number density, i.e., the 
number of particles per unit volume) [25], the bond- 
orientational order metric Qe the translational or¬ 
der metric T* [34, 39], the cumulative pore-size distri¬ 
bution F{S) [29], and the statistics of rattlers [39]. In 
addition, MRJ sphere packings exhibit disordered hy¬ 
peruniformity [29], meaning that they are locally disor¬ 
dered, but possess a hidden order on large length scales 
such that infinite-wavelength density fluctuations of MRJ 
packings vanish, i.e., the structure factor vanishes at 
the origin: \imk^oS{k) = 0 [43, 44]. Disordered hy¬ 
peruniformity can be seen as an “inverted critical phe¬ 
nomenon” with a direct correlation function c(r) that is 
long ranged [43, 45]. 

In this paper, we characterize the MRJ sphere pack¬ 
ings generated in Ref. [39] using Voronoi statistics, in¬ 
cluding certain types of correlation functions. We com¬ 
pare these computations to corresponding calculations 
for both a Poisson distribution of points and equilibrium 
hard-sphere liquids. In the second paper of this series, 
we will investigate density fluctuations, the pore-size dis¬ 
tribution, and two-point probability functions of MRJ 
packings. 

Many studies for disordered sphere packings have been 
devoted to computing the volume distribution of the 
Voronoi cells [e.g., 15, 17, 20, 46-52]; see Fig. 1 for a MRJ 
sphere packing and its Voronoi diagram [53]. However, 
such statistics are incomplete in that they only quantify 
local structural information. For example, with appro¬ 
priately rescaled variables, we will show that the distri¬ 
butions of the Voronoi volumes, surface areas, and inte¬ 
grated mean curvatures for the MRJ sphere packings are 
qualitatively similar to the distributions for an equilib¬ 
rium hard-sphere liquid and partly even for the spatially 
uncorrelated Poisson point process. 
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TABLE I. Mean (W^), standard deviation and correlation coefficients of the Minkowski functionals of single 

Voronoi cells in the Poisson point process, in a system of hard spheres in equilibrium at a packing fraction cf) = 0.48, and in the 
MRJ state. The unit of length is A = Ij , i.e., the number density p is set to unity. 




{W^) 


P/j.i 

Pm ,2 





Poisson 


Volume 

Wo 

1.0005(3) 

0.4230(2) 

0.98161(3) 

0.94486(8) 

Surface area 

Wi 

5.823(1) 

1.4798(7) 


0.98701(2) 

Integ. mean curv. 

W 2 

9.1623(8) 

1.0941(5) 






Equilibrium 


Volume 

Wo 

1.00000(7) 

0.07434(5) 

0.97700(5) 

0.9282(2) 

Surface area 

Wi 

5.4488(3) 

0.2681(2) 


0.98136(5) 

Integ. mean curv. 

W 2 

8.6477(2) 

0.2163(2) 

MRJ 


Volume 

Wo 

1.00000(3) 

0.04335(2) 

0.96976(4) 

0.9035(1) 

Surface area 

Wi 

5.4043(1) 

0.1695(8) 


0.97248(5) 

Integ. mean curv. 

W 2 

8.5894(1) 

0.1470(7) 




To quantify nonlocal structural information, we formu¬ 
late and compute correlation functions of the volume of 
Voronoi cells at a given distance and cell-cell probability 
density functions of finding a given sized Voronoi cell at 
a given distance of a sphere with another sized Voronoi 
cell. Because the volume is only one of a large class of 
versatile shape measures, namely, the Minkowski func¬ 
tionals [54-57], we devise and compute the correlation 
functions of all of the Minkowski functionals [58] . Besides 
characterizing MRJ packings in this way, we also carry 
out analogous calculations for the Poisson point process 
and the equilibrium hard-sphere liquid for purposes of 
comparison. We show that these Minkowski correlation 
functions contain visibly more information than the cor¬ 
responding standard pair-correlation functions, even in 
the case of the Poisson point process. 

In Sec. II, we analyze the distributions of the 
Minkowski functionals of the single Voronoi cells for 
the Poisson point process, equilibrium hard-sphere liq¬ 
uid configurations, and MRJ packings. In Sec. Ill, we 
define the aforementioned two different types of correla¬ 
tion functions. In Sec. IV, we determine the volume- 
volume correlation function numerically for the three- 
dimensional Poisson point process, equilibrium hard- 
sphere systems, and the MRJ state; we also calculate 
the correlation functions for the surface area, and the in¬ 
tegrated mean curvature. For a further investigation of 
the nonlocal structure features, we calculate in Sec. V the 
cell-cell probability density functions mentioned above. 
In Sec. VI, we make concluding remarks. 


II. MINKOWSKI FUNCTIONAL 
DISTRIBUTIONS OF A SINGLE VORONOI CELL 

While there are many detailed studies of the vol¬ 
ume distribution in disordered sphere packings [e.g., 
15, 17, 20, 46-52], here we analyze in a logarithmic plot 
the volume distributions of true MRJ packings as de¬ 
scribe above and extend the analysis to all three (non¬ 


trivial) Minkowski functionals: the volume, the surface 
area, and the integrated mean curvature. They are ro¬ 
bust and versatile shape descriptors which are widely 
used in statistical physics [55-57, 59] and in pattern anal¬ 
ysis [60-62]. We use vORO-|--t- [63, 64] to construct the 
Voronoi diagram of Poisson point patterns (about 1000 
patterns, each with 2000 points), equilibrium hard-sphere 
packings [3, 4] at a packing fraction (p = 0.48, which is 
slightly below the freezing transition (100 packings, each 
with 10000 spheres), and MRJ sphere packings produced 
by the TJ algorithm [25, 39] (about 1000 packings, each 
with 2000 spheres). The program karambola [65, 66] 
then computes the Minkowski functionals of each cell. 

We first determine the distributions of the three 
Minkowski functionals (Wq, volume; Wi, surface area; 
and W 2 , integrated mean curvature) of the single Voronoi 
cells. Table I provides the mean, the standard deviation, 
and the correlation coefficients Pau = 
of the three different Minkowski functionals. As a unit 
of length, we use A = Ij with p the number density, 
i.e., we compare the Poisson point process, the equilib¬ 
rium hard-sphere liquid, and the MRJ state at the same 
number density p = 1 (the unit volume contains one par¬ 
ticle on average) [67]. 

Because the number density is set to unity, the mean 
cell volume is also one. The average surface area and 
integrated mean curvature of a Voronoi cell in the MRJ 
state or in the equilibrium ensemble are slightly larger 
than those of a Poisson Voronoi cell because the latter 
is less regular, i.e., more aspherical. The Voronoi vol¬ 
ume fluctuations and the standard deviations of the other 
Minkowski functionals are much stronger in the irregular 
Poisson point process than in the hard-sphere packings, 
where the MRJ state has significantly smaller Voronoi 
volume fluctuations than the equilibrium hard-sphere liq¬ 
uid. The Minkowski functionals of a single Voronoi cell, 
e.g., its volume and its surface area, are strongly cor¬ 
related, i.e., the correlation coefficients p^^^ are at least 
0.9. The numerical estimates for the Poisson Voronoi tes- 
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FIG. 2. (Color online) Distributions of the single Minkowski 
functionals of a three-dimensional Voronoi cell in a Pois¬ 
son point process, an equilibrium hard-sphere system at a 
packing fraction (f) = 0.48, and a MRJ sphere packing: /r = 0 
volume (left), ^ = 1 surface area (center), and jj, — 2 
integrated mean curvature (right). The distributions are 
rescaled—like in Ref. [47] —by their mean (Ik)j) and their 
standard deviation criv^ (see Table I). The lines in the plot 
of the volume distributions are 7 distributions; generalized 7 
distributions are fitted to the distributions of the surface area 
and the integrated mean curvature. 


sellation are in agreement with the analytic values and 
numerical estimates in Ref. [ 68 ] and references therein. 

The high fidelity of the MRJ sphere packings pro¬ 
duced by the TJ algorithm allows one to study the rela¬ 
tion between the number of contacts of a sphere and the 
Minkowski functionals of its Voronoi cell. As expected, 
small cells have a higher number of contacts on average 
because a high local packing fraction [69] implies that 
there are many close neighbors. In units of A, the mean 
Voronoi volume of a rattler, i.e., an unjammed particle, is 
1.04 and that of a particle with 11 contacts is 0.88. The 
mean surface area of the Voronoi cells of rattlers and of 
backbone spheres with up to 11 contacts varies from 5.50 
to 4.92, respectively, and the average integrated mean 
curvature varies from 8.65 to 8.17, respectively. However, 
because of the small difference between near contacts and 
true contacts, the distributions of the Minkowski func¬ 
tionals for rattlers are only slightly shifted compared to 
the distributions of a typical cell. There are, for example, 
very small cells containing rattlers, which is consistent 
with previous findings [39]. 

Starr et al. [47] and, similarly, Aste et al. [15] showed 
that by shifting the volume distribution by its mean and 
rescaling with its standard deviation, the volume distri¬ 
butions of many different sphere packings collapse. Fig¬ 
ure 2 shows the rescaled distributions of the Minkowski 
functionals for the Poisson point process, the equilibrium 
hard-sphere liquid, and the MRJ packing. As expected, 
the volume distributions of the equilibrium hard-sphere 
packings and the MRJ packings are qualitatively very 
similar, while the distribution of the Poisson point pro¬ 
cess deviates. The same is true for the distribution of the 


mean curvatures. The distributions of the surface area for 
both the MRJ and the equilibrium hard-sphere packings 
are not only qualitatively similar to each other but also to 
the uncorrelated Poisson point process. So, besides the 
quantitative difference in the mean and the standard de¬ 
viations of the Minkowski functionals, the distributions 
of the Minkowski functionals of single Voronoi volumes 
are qualitatively similar for the equilibrium hard-sphere 
liquid and the MRJ state as well as even partially for the 
Poisson point process. The distribution of the Minkowski 
functionals of a single cell only incorporates local infor¬ 
mation and is rather insensitive to global structural fea¬ 
tures such as hyperuniformity of the MRJ state [29, 45]. 

Figure 3 shows the joint probability distribution of the 
volume and the surface area of a single Voronoi cell in 
a Poisson point process, an equilibrium hard-sphere liq¬ 
uid, and a MRJ sphere packing. The joint probability 
distributions for the equilibrium hard-sphere liquid and 
the MRJ state are also relatively similar. 

Both for the Poisson point process [70] and for 
many different numerical and experimental sphere pack¬ 
ings [15], the volume distribution follows well a 7 distri¬ 
bution [71]. We also find, for the volume distributions 
for the Poisson point process and the equilibrium hard- 
sphere liquid, an excellent agreement with 7 distribu¬ 
tions [72]. However, for the MRJ sphere packings there 
is a slight but statistically significant deviation for very 
small cells for which the frequency of occurrence is too 
high. The surface area and the integrated mean curvature 
distributions are well approximated by generalized 7 dis¬ 
tributions [73], which was already found for the Poisson 
point process by Refs. [51, 74]. However, the distribu- 
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FIG. 3. (Color online) Joint distribution of volume v and sur¬ 
face area S' of a single three-dimensional Voronoi cell in a log¬ 
arithmic scale for a Poisson point process and a MRJ sphere 
packing; the joint distribution for the equilibrium hard-sphere 
liquid is represented by the blue contour plot. The unit of 
length is A = 1/p^f^, where p is the number density. Samples 
of a MRJ sphere packing and an overlapping sphere packing, 
where the sphere centers follow a Poisson point process, are 
depicted together with their Voronoi diagrams. 
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tions for the MRJ sphere packings deviate slightly but 
statistically significantly from a generalized 7 distribu¬ 
tion for cells with small surface area or small integrated 
mean curvature, respectively [75]. 


III. CORRELATION FUNCTIONS AND 
PROBABILITY DENSITY FUNCTIONS OF 
MINKOWSKI FUNCTIONALS 


In order to quantify the global structure of the Voronoi 
diagram, correlation functions of the Minkowski function¬ 
als of cells at a distance r and cell-cell probability density 
functions are introduced and defined here. 


A. Correlation Functions of Minkowski Functionals 

We define the volume-volume correlation function 
Coo(’'i) *’ 2 ) of the Voronoi cells of an arbitrary point pro¬ 
cess as the correlation between the volume of two Voronoi 
cells given that the corresponding centers are at the po¬ 
sitions ri and r 2 '. 

^ ... N {viri)v{r2)) - {v{ri)){v(r2)) ,,, 

Coo{ri,r2) ■■= — -———( 1 ) 

'^v{ri\r2)^v{r2\ri) 

where (.) denotes the ensemble average given two points 
at Ti and r 2 ] and is the standard deviation of 

the volume v of the Voronoi cell at ri given that there is 
another point at rj. Note that because of this condition, 
both the mean and the standard deviation of a single 
Voronoi volume are functions of the positions ri and r 2 '. 
e.g., knowing that there is a point in close proximity, 
very large volumes are less likely and the mean volume 
decreases. For a statistically homogeneous and isotropic 
point process, the volume-volume correlation is simply 
a radial function, which we denote by C'oo(c), where 
1^2 —J'lli- The correlation function Coo(r) S [—1;1] 
measures the correlations, both positive and negative 
(anticorrelations), between Voronoi volumes of cells given 
that their centers are at a distance r. 

The Voronoi tessellation assigns to each point a volume 
of its corresponding Voronoi cell. This is a special case 
of a marked point process where the constructed mark 
assigned to each point is determined by the positions of 
the points in the neighborhood. In this sense, the volume- 
volume correlation function can be seen as a special type 
of a marked correlation function [56, 68 , 76]. 

The volume-volume correlation function does not, in 
general, converge to perfect correlation for vanishing ra¬ 
dial distance limr_>.o C'oo (’’) < 1 because for all r > 0 the 
correlation function Cqo (r) provides the correlation of the 
Voronoi volumes of two different cells with volumes u(0) 
and v{r). Because the cell is perfectly correlated with 
itself, i.e., C'oo(O) = 1, the correlation function Coo(r) 
is discontinuous at the origin. If there is no long-range 
order, the correlation function tends to zero for infinite 
radial distance limj._>oo C'oo(^) — 0 . 


The correlation functions of the other Minkowski func¬ 
tionals are defined analogously to Eq. (1), replacing vol¬ 
ume (/J, = 0 ) by surface area (p, = 1 ) or integrated mean 
curvature (^ = 2 ): 

^ ^ {Wi2iri)W^{r2))-{W^{ri)){W^{r2)) 

(r 1 1 r 2 ) (r2 1 r 1) 

( 2 ) 

with a\Y^(^ri\rj) the standard deviation of the Minkowski 
functional IV^ of the Voronoi cell at given that there is 
another point at rj. For a statistically homogeneous and 
isotropic point process, the correlation function of the 
Minkowski functionals is again a radial function, which 
we denote by C^^{r). In general, will be discon¬ 

tinuous for r — 0 , as noted above for the volume-volume 
correlation function. 

In the Appendix, we calculate the volume-volume cor¬ 
relation function analytically for the one-dimensional 
Poisson point process. In Sec. IV, we determine the 
correlation functions for the three-dimensional Poisson 
point process, the equilibrium hard-sphere liquid, and 
MRJ sphere packings. 

A different type of correlation function, a pointwise 
Voronoi correlation function, assigns to arbitrary points 
the volumes of the Voronoi cells in which they lie [52]. 
Correlations between Voronoi volumes have also already 
been studied by finding a nonlinear scaling in the ag¬ 
gregate Voronoi volume fluctuations as a function of the 
sample size [49]. 


B. Cell-Cell Probability Density Functions of the 
Voronoi Volume 

The volume-volume correlation function Coo{r 2 ,ri) is 
defined conditionally on the fact that the centers of the 
two cells are at Ti and r 2 . The full two-point informa¬ 
tion about the Voronoi volumes is given by the cell-cell 
probability density function p(t’ 2 , u, ri, v*) of finding two 
points in the point process at two arbitrary positions r 2 
and Vi with associated Voronoi volumes v and u*, respec¬ 
tively. It quantifies, for example, how likely it is to find 
near a point with a small Voronoi cell another point with 
either a large or another small Voronoi cell. Integrat¬ 
ing over the volumes yields the standard pair-correlation 
function. 


52 (^ 2 , Ti) = dvdv 


p(r2,u,ri,u*) 

p(r- 2 )p(r'i) 


( 3 ) 


This relation clearly indicates that the Minkowski prob¬ 
ability density function p{r 2 ,v,ri,v*) contains more in¬ 
formation than Moreover, the volume-volume 

correlation function C'oo('i' 2 , »’i) from Sec. Ill A follows 
from calculating the moments (vv*), (v), and (v*) of 
p{r 2 ,v,ri,v ) —^ corresponding standard devi- 

p{r2)p{ri)g2{r2,ri) ^ » 

ations ay and ay*. 
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For a statistically homogeneous and isotropic point 
process, the cell-cell probability density function is a ra¬ 
dial function, i.e., it only depends on the radial distance 
r = |t '2 — ri |: p(r, v,v*) is the probability density of find¬ 
ing two points with Voronoi volumes v and v* at a radial 
distance r. If there is no long-range order, the cell-cell 
probability density function p{r, v, v*) converges for large 
radii r —> oo to f{v)f{v*), with p the number density 
and f{v) the distribution of the Voronoi volume u of a 
single cell (see Sec. II). 

For a better visualization and comparison of different 
volumes, we divide the cell-cell probability density func¬ 
tion by its long-range value; the cell-cell pair-correlation 
function is defined as 


gvv{r2,v,ri,v*) 


p{r2,v,ri,v*) 

p{r2)fiv)piri)f{v*) 


( 4 ) 


and, for a homogeneous and isotropic system. 


gvv{r,v,v*) 


Pir,v,v*) 
P^fiv)fiv*) ■ 


( 5 ) 


If gyy{r,v,v*) > I, it is more likely to find a pair of 
Voronoi cells with volumes v and u* at a distance r than 
to find them at a large distance, i.e., uncorrelated. If 
gvvir,v,v*) < I, the occurrence of a point in the point 
process with a Voronoi volume u at a distance r of an¬ 
other Voronoi center with a Voronoi volume v* is sup¬ 
pressed. Analogous cell-cell pair-correlation functions 
can be defined for the other Minkowski functionals. 

We analytically calculate the cell-cell probability den¬ 
sity function for the one-dimensional Poisson point pro¬ 
cess in the Appendix. In Sec. V, we determine the cell-cell 
pair correlation function for the three-dimensional Pois¬ 
son point process, the equilibrium hard-sphere liquid, and 
MRJ sphere packings. 


IV. CORRELATION FUNCTIONS OF 
MINKOWSKI FUNCTIONALS 


In order to sample the correlation functions of the 
Minkowski functionals, the distances of all pairs of par¬ 
ticles [77] are computed and assigned to a bin. For each 
radial distance, i.e., for each bin, the correlation coef¬ 
ficient of the Minkowski functionals of the two Voronoi 
cells is determined. 

Figures 4-6 compare the correlation functions of the 
Minkowski functionals for the Poisson point process, 
equilibrium hard-sphere liquids, and MRJ sphere pack¬ 
ings. It is seen that these Minkowski correlation functions 
contain visibly more information than the corresponding 
standard pair-correlation functions, even in the case of 
the Poisson point process. 



FIG. 4. (Color online) Correlation functions for a three- 
dimensional Poisson point process', pair-correlation function 
g 2 {r) (bottom); volume-volume correlation function C'oo(r) 
(center), which is the correlation function of the volumes 
of two Voronoi cells given that their centers are at a dis¬ 
tance r; the mark correlation functions of the three different 
Minkowski functionals (top): p = 0 volume, p = 1 surface 
area, and p = 2 integrated mean curvature. The radial dis¬ 
tance r is normalized by A = Ijp^^^, where p is the number 
density. 


A. Poisson Point Process 


It is evident that in the infinite-system limit, the pair- 
correlation function g 2 {r) is a constant (unity) for the 
Poisson point process, i.e., the points are completely un¬ 
correlated. Because a Voronoi cell is determined by the 
neighbors of its center, the volume will obviously be cor¬ 
related; see Fig. 4. There are large Voronoi volume fluc¬ 
tuations for the Poisson point process. Very large cells 
lead to a strong correlation of the Voronoi volumes even 
for distances up to four times the mean nearest-neighbor 
distance. This is to be contrasted with the standard pair- 
correlation function g 2 if), which is trivially unity for all 
radial distances. 

Figure 4 compares the correlation functions for 

all Minkowski functionals p — 0,1, 2. All functionals have 
approximately the same correlation length. For r —>■ 0, 
the surface areas are more strongly correlated than the 
volumes because at small radial distances, the cells will 
most likely share a face. In the Appendix, we calculate 
C'oo(f) analytically for the one-dimensional Poisson point 
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FIG. 5. (Color online) Correlation functions for equilibrium 
hard spheres with diameter D at a packing fraction 0 = 0.48; 
for details, see Fig. 4. 


process. 

B. Equilibrium Hard-Sphere Liquid 

Figure 5 shows the pair-correlation function and the 
correlation functions of the Minkowski functionals for 
equilibrium hard-sphere liquid configurations at a pack¬ 
ing fraction (j) = 0.48. Because the hard spheres are im¬ 
penetrable, the correlation functions of the Minkowski 
functionals are only defined for radial distances larger 
or equal to the diameter D of a sphere; in this case, 
D « 0.97 A. 

There is a strong correlation of the Voronoi volumes 
of spheres that are in near contact because the Voronoi 
neighbors are correlated by construction of the Voronoi 
diagram. However, the maximum correlation is reached 
for noncontacting spheres at r « 1.3 A; a large cell has 
many neighbors and a Voronoi neighbor with a sphere 
not in contact will be, on average, larger than another 
neighbor cell with a contacting sphere. 

Between 1.8 A and 2.4 A, there is a double peak of an¬ 
ticorrelation and, for larger radial distances, there is an 
oscillating anticorrelation and correlation similar to the 
pair correlation function ^ 2 , but nearly inverted. The 
correlation length of the Voronoi volumes in the hard- 
sphere liquid is larger than in the uncorrelated Poisson 
point process, where the correlation was only due to the 



FIG. 6 . (Color online) Correlation functions for MRJ sphere 
packings of spheres with diameter D; the average packing frac¬ 
tion is (j) ~ 0.64. For details, see Fig. 4. The pair-correlation 
function 32 (r) is in agreement with previous results for MRJ 
sphere packings [39]. 

large Voronoi volume fluctuations. 

At the top of Fig. 5, the correlation functions of the 
other Minkowski functionals are compared. Similar to 
the Poisson case, the integrated mean curvature is more 
strongly correlated at contact r = D than the surface 
area which, in turn, is more strongly correlated than the 
volume. There is no double anticorrelation peak in the 
integrated mean curvature. For large radii, the correla¬ 
tion functions are shifted against each other despite the 
strong correlation of the different functionals for a single 
Voronoi cell. The surface area and the integrated mean 
curvature are slightly less (anti)correlated. 

C. MRJ Sphere Packings 

The pair-correlation function g 2 (r) and the correlation 
functions of the Minkowski functionals of the MRJ sphere 
packings are shown in Fig. 6. The diameter of the spheres 
in the MRJ sphere packings is ZJ « 1.07A. The most 
striking differences in the pair correlation of the jammed 
packings to the equilibrium packings are the two discon¬ 
tinuities at r = \f^D and r = 2D, the split-second peak, 
which corresponds to configurations of two edge-sharing 
equilateral and coplanar triangles (r = VSD) or a linear 
chain of three particles (r = 2D ), respectively [78] . There 






























is also a significant (seemingly nonanalytic) feature of the 
volume-volume correlation function Coo(r) at r = \/3D: 
a dip in the anticorrelation. However, at r = 2D, the 
feature is statistically insignificant. At least two double 
anticorrelation peaks are clearly resolved. 

The most important qualitative differences in the 
volume-volume correlation function are the much 
stronger anticorrelations in the MRJ packings compared 
to the equilibrium packings. The correlation with the 
nearest neighbors is weaker and the first anticorrelation 
double peak is more than twice as strong as for the equi¬ 
librium hard-sphere packings. The MRJ sphere pack¬ 
ings are hyperuniform [29, 45], i.e., large-scale density 
fluctuations are suppressed. Therefore, strong Voronoi 
volume anticorrelations are necessary such that Voronoi 
cells with a high local packing fraction are accompanied 
by cells with rather low packing fractions, and vice versa. 

Another difference between MRJ and equilibrium 
packings is a stronger shift of the correlation func¬ 
tions of the other Minkowski functionals. For the 
MRJ packings, there are radial distances, e.g., r = 
2.51 A, at which the integrated mean curvatures are anti¬ 
correlated [(722(2.51 A) < 0] but the volumes are corre¬ 
lated [(7oo(2.51 A) > 0], and vice versa. 

So, in contrast to the local Voronoi analysis, the global 
Voronoi analysis of the MRJ packing reveals qualita¬ 
tive structural differences to the equilibrium hard-sphere 
liquid. 


V. CELL-CELL PROBABILITY DENSITY 
FUNCTIONS 

The sampling of the cell-cell probability density func¬ 
tion p{r,v,v*) is very similar to that of the pair- 
correlation function 52 (see, e.g.. Ref. [4]); only an ad¬ 
ditional binning with respect to the Voronoi volumes is 
needed. Figures 7-9 show the cell-cell pair-correlation 
function = p{r,v,v*)/(p^f{v)f(v*)) for exemplary 
large or small cell volumes v,v* in the three-dimensional 
Poisson point process, in an ensemble of equilibrium hard 
spheres, or in the MRJ sphere packings. As examples of 
large or small cells, the volumes were chosen such that 
their probability density is equal to 1 /3 of the maximum 
of the volume distribution; see Table I and Fig. 2. 

For the Poisson point process, the small cells are 
strongly correlated at short distances because, by con¬ 
struction, there must be points at close distances, and 
the neighbor cells of a small Voronoi cell are more likely 
to be small as well. However, the probability of finding 
a point with either a corresponding large or a small cell 
at a short radial distance of the center of a large cell is 
strongly suppressed because it is unlikely for the center 
of a large cell to have close neighbors. At intermediate 
distances, two large cells are correlated, as expected, be¬ 
cause of the Voronoi construction. 

In the equilibrium hard-sphere liquid, the large cells 
at near contact are less correlated than the small cells. 



FIG. 7. (Color online) The cell-cell pair-correlation function 
gvv{r,v,v*) for a three-dimensional Poisson point process for 
either two large cells (bottom), a large and a small cell (cen¬ 
ter), or two small cells (top). As examples of large or small 
cells, the volumes were chosen such that their probability den¬ 
sity is equal to 1/3 of the maximum of the volume distribution; 
see Table I and Fig. 2. The curves are compared to the stan¬ 
dard pair-correlation function 92 (r) (dashed black line), which 
is trivially unity for the uncorrelated Poisson point process. 
The radial distance r is normalized by A = where p is 

the number density. 

However, at slightly larger distances, where 52 shows an¬ 
ticorrelation and the small cells are even more strongly 
anti-correlated, the large cells are positively correlated. 
For distances larger than twice the diameter, the cell¬ 
cell pair-correlation function for a large and a small cell 
is equal to the standard pair-correlation function within 
statistical significance. However, both the cell-cell pair- 
correlation functions of finding two short or of finding two 
large cells at large radial distance r are shifted compared 
to the standard pair-correlation function. 

These features can also be found in the MRJ sphere 
packings. Moreover, the anticorrelations of two small 
cells are much stronger. The split-second peak even 
separates in two stronger peaks with anticorrelation in 
between, where the standard pair-correlation function 
shows positive correlation. In contrast to this, the peak 
at r = completely vanishes for two large cells 

<7j;«('\/ 37A, 1.06,1.06) = 1, and the peak at r = 2D is 
significantly weaker. 


VI. CONCLUSIONS AND DISCUSSIONS 

We have characterized the structure of MRJ sphere 
packings by computing the Minkowski functionals, i.e.. 
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FIG. 8. (Color online) The cell-cell pair-correlation function 
gvv{r,v,v*) for equilibrium hard spheres at a global packing 
fraction (j) = 0.48; for details, see Fig. 7. 


FIG. 9. (Color online) The cell-cell pair-correlation function 
gvv{r,v,v*) of the MRJ sphere packings; the average packing 
fraction is (j> ks 0.64; for details, see Fig. 7. 


the volume, the surface area, and the integrated mean 
curvature, of the associated Voronoi cells. The local anal¬ 
ysis, i.e., the probability distribution of the Minkowski 
functionals of a single Voronoi cell, provides qualitatively 
similar results for the equilibrium hard-sphere liquid and 
the MRJ packings and partly even for the uncorrelated 
Poisson point process. 

In order to study the global structure of the Voronoi 
cells, we have improved upon this analysis by introduc¬ 
ing the correlation functions of the Minkowski 

functionals and the cell-cell probability density func¬ 
tion p{r,v,v*). The correlation function C^^{r) mea¬ 
sures the correlation of the Minkowski functionals 
of two Voronoi cells given that the corresponding cen¬ 
ters are at a distance r. The cell-cell probability den¬ 
sity function p(r,v,v*) also incorporates the probability 
that there are two particles at a distance r. For an eas¬ 
ier interpretation and better visualization, we have de¬ 
fined the dimensionless cell-cell pair-correlation function 
gw = p{r,v,v*)/{p'^f{v)f{v*)), where f{v) is the proba¬ 
bility of the Voronoi volume v. The generalization of the 
pair-correlation to the cell-cell pair correlations provides 
powerful theoretical and computational tools to charac¬ 
terize the complex local geometries that arise in jammed 
disordered sphere packings. 

Because the faces of a Voronoi cell are bisections be¬ 
tween a point in the point process (whether a packing or 
not) and its neighbors and, moreover, because Voronoi 
neighbors share a face and edges, the Minkowski func¬ 
tionals of neighboring Voronoi cells are correlated by 
construction. This leads to a large correlation length 
for the Voronoi cells in a Poisson point process because 
of large Voronoi volnme fluctuations. In the equilib¬ 


rium hard-sphere liquid and MRJ sphere packings, there 
are correlations and anticorrelations. In contrast to the 
qualitatively similar local Voronoi structure, the global 
Voronoi structure of the MRJ hard-sphere packings is 
qualitatively quite different from that of an equilibrium 
hard-sphere liquid. We find strong Voronoi volume anti¬ 
correlations, which is consistent with previous findings 
that MRJ sphere packings are hyperuniform [29, 45], 
i.e., large-scale density flnctuations are suppressed. MRJ 
sphere packings are prototypical glasses in that they have 
no long-range order but they are perfectly rigid, i.e., the 
elastic moduli are unbounded [19, 37, 79]. The global 
analysis introduced here shows the difference in the struc¬ 
ture of the Voronoi cells of the MRJ state and those of a 
hard-sphere liquid, which further indicates that the struc¬ 
ture of a glass is not that of a “frozen liquid” [45, 79, 80]. 

An already known distinct strnctural difference be¬ 
tween the hyperuniform MRJ sphere packings and equi¬ 
librium hard-sphere liquids is that while in the equilib¬ 
rium packing the total pair-correlation function h{r) = 
52 ("c) — 1 is exponentially damped, the total correlation 
function of the MRJ state has a negative algebraic power- 
law tail [29, 45, 79]. It is an interesting question as to 
whether the asymptotic behaviors of the correlation func¬ 
tion of the Minkowski functionals C'^^(r) or the radial 
cell-cell correlation functions Pwir, v, v*) are different for 
the MRJ state and the hard-sphere liquid. However, a 
direct observation of the power-law tail has, so far, not 
been possible [29, 45, 79]; much larger systems are needed 
but are not available at the moment. Still, the global 
characteristics C'^^(r) and 5 „„(r, v,v*), introduced in the 
present paper, allow for an investigation of the underlying 
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geometrical reasons for the negative algebraic tail in the 
total pair-correlation function: the suppressed clustering 
of regions with low and high local packing fractions [45]. 

Moreover, they also allow for a quantification of the 
global structure of other cellular structures, e.g., foams, 
where the centers of mass of the single cells can be used 
as centers of the cells instead of the Voronoi centers used 
here. 

A frequently discussed question is whether or not there 
are local icosahedral configurations in jammed pack¬ 
ings [14, 20, 21, 31, 81], i.e., a central sphere with 12 
spheres in contact where the centers of the touching 
spheres form a regular icosahedron. The Voronoi cell 
of the central sphere in such an icosahedron is a regu¬ 
lar dodecahedron, which has the maximum possible local 
packing fraction (« 0.76). There is growing evidence, 
that there are no regular icosahedral arrangements in 
hard-sphere packings, e.g., see Refs. [14, 78, 82]. In¬ 
deed, we find in our MRJ sphere packings no regular and 
hardly any nearly regular dodecahedral Voronoi cells. All 
spheres out of more than two million have less than 12 
contacts. There are local packing fractions up to 0.75, 
but only 4.2 x 10“® of all cells have a local packing frac¬ 
tion greater than 0.74. 

In a preliminary approach to look for possibly strongly 
distorted dodecahedral Voronoi cells in the MRJ sphere 
packings, we examined the topology of the Voronoi poly- 
hedra, i.e., the number of faces and the corresponding 
types of polygons, following Refs. [81, 83, 84]. In a com¬ 
pact notation, the topology of a polyhedron is given by 
the so-called p vector (77,3 714 715 ng), where 713 is the num¬ 
ber of triangles, 774 of quadrilaterals, 775 of pentagons, 
and Tig of hexagons. The dodecahedron is formed by 12 
pentagons, i.e., its topology is denoted by (0 0 12 0). Al¬ 
though these polyhedron characteristics are discontinu¬ 
ous and inadequately metric for definite conclusions [85] , 
they can provide a first insight into whether there could 
be a significant number of distorted dodecahedra. In the 
MRJ sphere packings, 1.1% of all cells have the topol¬ 
ogy of a dodecahedron (0 0 12 0) [86]. The average local 
packing fraction of those distorted dodecahedra is 0.69 
and is thus significantly greater than the total mean lo¬ 
cal packing fraction which is 0.64. However, only 0.4 % of 
the distorted dodecahedra have a local packing fraction 
greater than 0.74. The distorted dodecahedra also have a 
higher average number of contacts, « 7, compared to the 
typical cell, « 6, but as stated above there is not a single 
sphere with 12 contacts in this high-quality MRJ data. 
There are 25 other topologies in the Voronoi diagram 
of the MRJ sphere packings that occur more frequently 
than the dodecahedron. With 5.2 % of all cells, the most 
likely topology is (0 3 6 5). However, by adding one or 
two faces, the dodecahedron can transform to the follow¬ 
ing polyhedra [81]: 1.1% of all cells in the MRJ sphere 
packings are (1 0 9 3), 3.1 % are (0 1 10 2), and 4.4% are 
(0 2 8 4). The latter is the second most common type in 
the MRJ sphere packings. So, while we find no regular 
icosahedral configurations in the MRJ sphere packings. 


the preliminary topological analysis indicates that more 
detailed studies of probably strongly distorted icosahe- 
dra could be interesting. For example, also in metallic 
glasses, significant numbers of distorted icosahedra have 
been found [87, 88]. 

In the second paper of this series, we will further inves¬ 
tigate the global structure of the MRJ sphere packings by 
looking at density fluctuations, the pore-size distribution, 
and the two-point correlation functions. 
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Appendix A: Correlation Fnnctions of the Voronoi 
cells in the One-Dimensional Poisson Point Process 


As an introductory example of the global Voronoi 
statistics introduced in this paper, we analytically calcu¬ 
late the correlation functions of the Voronoi cells in the 
one-dimensional Poisson point process. Therefore, we use 
the probability density functions Hp(ni) and Hp{nr) of 
the nearest neighbor on the left-hand side at a distance 
71; or on the right-hand side at a distance Ur, respectively. 

Two points xi and X 2 at a distance r are given. 
Without loss of generality, we assume in the following 
xi = Q and X 2 = r. The nearest-neighbor probabil¬ 
ity density functions of xi are Hp(ni) = pe“^"' and 
Hp{nr) = pe~P'^''9{r — Ur) + e~P'^5(r — Ur), where 6{x) is 
the Heaviside step function and 6{x) is the Dirac delta 
distribution. For X 2 , the distributions for the right- and 
the left-hand side simply exchange. The probability dis¬ 
tribution f{v*\r) of the volume v* of the cell correspond¬ 
ing to Xi = 0 is given by the average of i5( "^+"i — v*): 


Given a volume v* of the cell corresponding to Xi: If 
V* < rl2, there will be at least one additional point y 
between xi and X 2 - Its distance z to X 2 is uniformly dis¬ 
tributed between r — 2v* and r. With h(z\r,v*) denoting 
the probability density function of this distance, the con¬ 
ditional probability distribution of the volume v of the 
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FIG. 10. (Color online) Correlation functions for the one¬ 
dimensional Poisson point process: pair-correlation function 
g 2 {r) (bottom); volume-volume correlation function Coo{i~) 
(top). The distance r is scaled (in this one-dimensional ex¬ 
ample) by the inverse of the number density A = 1/p. An 
example of a one-dimensional Voronoi diagram of a Poisson 
point process is depicted. 


cell corresponding to X 2 is given by 

f[v\r,v'')= [ dzhizlr.v'^) • f{v\z) (A2) 

with f{v\z) from Eq. (Al). A case-by-case analysis for 
differently large v compared to r and v* is needed. If 

f{v\r,v*)=4vp'^e~^'"^. (A3) 

E^-v* <v<^, then 

f{v\r, V*) = [4:v{r - v)p - (r - 2v*fp 

-|-4(w -b V*) — 2r] . 

If i; > then 

fiv\r, V*) = 2pe-^'’P [1 + p{r - ?;*)] . (A5) 

If ri* > r/2, there will be at least one additional point y 
between xi and X 2 with probability rp/(rp+ 1) and with 
probability l/(rp-|- 1) the points Xi and X 2 are nearest 


neighbors. In the first case, the conditional probability 
distribution of the volume v is given by Eq. (A2), but the 
distance z is now uniformly distributed between 0 and r. 
If u < then 


4une 

= ^ - ((r-v)p+l). (A6) 

r 

If V > then 


f{v\r,v*) = pe '^^P{rp + 2). (A7) 

In the second case, where there is no point between Xi 
and X 2 ^ the volume v is at least r/2 and completely de¬ 
termined by the nearest neighbor of X 2 on the right-hand 
side, 


f{v\r,v*) = 2pe''Pe-'^'’P. (A8) 

The cell-cell probability density function p(r,v,v*) 
from Sec. Ill B is then given by 

p{r,v,v*) = p^ ■ f{v*\r) ■ f[v\r,v *). (A9) 

Erom Eqs. (Al) and (A3) follows the asymptotic behavior 
of p(r, V, V*) > p'^f{v)f(v*). 

As described in Sec. IIIB, the volume-volume corre¬ 
lation function Coo from Sec. Ill A follows straightfor¬ 
wardly. 


Cm{r) 


r^p^ -2rp + 2- 2e-^P 
4 - 4rpe-’'^' - 2e-2’'P 


(AlO) 


Figure 10 depicts the volume-volume correlation func¬ 
tion Coo(r) of the one-dimensional Poisson point process; 
both the analytic result and simulation data are shown. 

As discussed in Sec. IV A for the three-dimensional 
Poisson point process, the Voronoi neighbors are corre¬ 
lated by construction. Although very large Voronoi cells 
are rather unlikely, their next-neighbor correlation leads 
to a large correlation length in Coo(r). In contrast to the 
three-dimensional case, the Voronoi neighbors are uncor¬ 
related if the distance of their centers vanishes because in 
one dimension these Voronoi cells become independent. 
They only depend on either the nearest neighbor on the 
left- or on the right-hand side of xi = X 2 , which are in¬ 
dependent of each other. For large radii, the correlation 
vanishes exponentially, as expected, because there is no 
long-range order in the Voronoi diagram. 
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